home *** CD-ROM | disk | FTP | other *** search
/ The Atari Compendium / The Atari Compendium (Toad Computers) (1994).iso / files / prgtools / editors / emcs1857 / 1857el~1.zoo / lisp / float.el < prev    next >
Encoding:
Text File  |  1992-02-03  |  14.7 KB  |  452 lines
  1. ;; Copyright (C) 1986 Free Software Foundation, Inc.
  2. ;; Author Bill Rosenblatt
  3.  
  4. ;; This file is part of GNU Emacs.
  5.  
  6. ;; GNU Emacs is free software; you can redistribute it and/or modify
  7. ;; it under the terms of the GNU General Public License as published by
  8. ;; the Free Software Foundation; either version 1, or (at your option)
  9. ;; any later version.
  10.  
  11. ;; GNU Emacs is distributed in the hope that it will be useful,
  12. ;; but WITHOUT ANY WARRANTY; without even the implied warranty of
  13. ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14. ;; GNU General Public License for more details.
  15.  
  16. ;; You should have received a copy of the GNU General Public License
  17. ;; along with GNU Emacs; see the file COPYING.  If not, write to
  18. ;; the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  19.  
  20. ;; Floating point arithmetic package.
  21. ;;
  22. ;; Floating point numbers are represented by dot-pairs (mant . exp)
  23. ;; where mant is the 24-bit signed integral mantissa and exp is the
  24. ;; base 2 exponent.
  25. ;;
  26. ;; Emacs LISP supports a 24-bit signed integer data type, which has a
  27. ;; range of -(2**23) to +(2**23)-1, or -8388608 to 8388607 decimal.
  28. ;; This gives six significant decimal digit accuracy.  Exponents can
  29. ;; be anything in the range -(2**23) to +(2**23)-1.
  30. ;;
  31. ;; User interface:
  32. ;; function f converts from integer to floating point
  33. ;; function string-to-float converts from string to floating point
  34. ;; function fint converts a floating point to integer (with truncation)
  35. ;; function float-to-string converts from floating point to string
  36. ;;                   
  37. ;; Caveats:
  38. ;; -  Exponents outside of the range of +/-100 or so will cause certain 
  39. ;;    functions (especially conversion routines) to take forever.
  40. ;; -  Very little checking is done for fixed point overflow/underflow.
  41. ;; -  No checking is done for over/underflow of the exponent
  42. ;;    (hardly necessary when exponent can be 2**23).
  43. ;; 
  44. ;;
  45. ;; Bill Rosenblatt
  46. ;; June 20, 1986
  47. ;;
  48.  
  49. (provide 'float)
  50.  
  51. ;; fundamental implementation constants
  52. (defconst exp-base 2
  53.   "Base of exponent in this floating point representation.")
  54.  
  55. (defconst mantissa-bits 24
  56.   "Number of significant bits in this floating point representation.")
  57.  
  58. (defconst decimal-digits 6
  59.   "Number of decimal digits expected to be accurate.")
  60.  
  61. (defconst expt-digits 2
  62.   "Maximum permitted digits in a scientific notation exponent.")
  63.  
  64. ;; other constants
  65. (defconst maxbit (1- mantissa-bits)
  66.   "Number of highest bit")
  67.  
  68. (defconst mantissa-maxval (1- (ash 1 maxbit))
  69.   "Maximum permissable value of mantissa")
  70.  
  71. (defconst mantissa-minval (ash 1 maxbit)
  72.   "Minimum permissable value of mantissa")
  73.  
  74. (defconst floating-point-regexp
  75.   "^[ \t]*\\(-?\\)\\([0-9]*\\)\
  76. \\(\\.\\([0-9]*\\)\\|\\)\
  77. \\(\\(\\([Ee]\\)\\(-?\\)\\([0-9][0-9]*\\)\\)\\|\\)[ \t]*$"
  78.   "Regular expression to match floating point numbers.  Extract matches:
  79. 1 - minus sign
  80. 2 - integer part
  81. 4 - fractional part
  82. 8 - minus sign for power of ten
  83. 9 - power of ten
  84. ")
  85.  
  86. (defconst high-bit-mask (ash 1 maxbit)
  87.   "Masks all bits except the high-order (sign) bit.")
  88.  
  89. (defconst second-bit-mask (ash 1 (1- maxbit))
  90.   "Masks all bits except the highest-order magnitude bit")
  91.  
  92. ;; various useful floating point constants
  93. (setq _f0 '(0 . 1))
  94.  
  95. (setq _f1/2 '(4194304 . -23))
  96.  
  97. (setq _f1 '(4194304 . -22))
  98.  
  99. (setq _f10 '(5242880 . -19))
  100.  
  101. ;; support for decimal conversion routines
  102. (setq powers-of-10 (make-vector (1+ decimal-digits) _f1))
  103. (aset powers-of-10 1 _f10)
  104. (aset powers-of-10 2 '(6553600 . -16))
  105. (aset powers-of-10 3 '(8192000 . -13))
  106. (aset powers-of-10 4 '(5120000 . -9))
  107. (aset powers-of-10 5 '(6400000 . -6))
  108. (aset powers-of-10 6 '(8000000 . -3))
  109.  
  110. (setq all-decimal-digs-minval (aref powers-of-10 (1- decimal-digits))
  111.       highest-power-of-10 (aref powers-of-10 decimal-digits))
  112.  
  113. (defun fashl (fnum)            ; floating-point arithmetic shift left
  114.   (cons (ash (car fnum) 1) (1- (cdr fnum))))
  115.  
  116. (defun fashr (fnum)            ; floating point arithmetic shift right
  117.   (cons (ash (car fnum) -1) (1+ (cdr fnum))))
  118.  
  119. (defun normalize (fnum)
  120.   (if (> (car fnum) 0)            ; make sure next-to-highest bit is set
  121.       (while (zerop (logand (car fnum) second-bit-mask))
  122.     (setq fnum (fashl fnum)))
  123.     (if (< (car fnum) 0)        ; make sure highest bit is set
  124.     (while (zerop (logand (car fnum) high-bit-mask))
  125.       (setq fnum (fashl fnum)))
  126.       (setq fnum _f0)))            ; "standard 0"
  127.   fnum)
  128.       
  129. (defun abs (n)                ; integer absolute value
  130.   (if (natnump n) n (- n)))
  131.  
  132. (defun fabs (fnum)            ; re-normalize after taking abs value
  133.   (normalize (cons (abs (car fnum)) (cdr fnum))))
  134.  
  135. (defun xor (a b)            ; logical exclusive or
  136.   (and (or a b) (not (and a b))))
  137.  
  138. (defun same-sign (a b)            ; two f-p numbers have same sign?
  139.   (not (xor (natnump (car a)) (natnump (car b)))))
  140.  
  141. (defun extract-match (str i)        ; used after string-match
  142.   (condition-case ()
  143.       (substring str (match-beginning i) (match-end i))
  144.     (error "")))
  145.  
  146. ;; support for the multiplication function
  147. (setq halfword-bits (/ mantissa-bits 2)    ; bits in a halfword
  148.       masklo (1- (ash 1 halfword-bits)) ; isolate the lower halfword
  149.       maskhi (lognot masklo)        ; isolate the upper halfword
  150.       round-limit (ash 1 (/ halfword-bits 2)))
  151.  
  152. (defun hihalf (n)            ; return high halfword, shifted down
  153.   (ash (logand n maskhi) (- halfword-bits)))
  154.  
  155. (defun lohalf (n)            ; return low halfword
  156.   (logand n masklo))
  157.  
  158. ;; Visible functions
  159.  
  160. ;; Arithmetic functions
  161. (defun f+ (a1 a2)
  162.   "Returns the sum of two floating point numbers."
  163.   (let ((f1 (fmax a1 a2))
  164.     (f2 (fmin a1 a2)))
  165.     (if (same-sign a1 a2)
  166.     (setq f1 (fashr f1)        ; shift right to avoid overflow
  167.           f2 (fashr f2)))
  168.     (normalize
  169.      (cons (+ (car f1) (ash (car f2) (- (cdr f2) (cdr f1))))
  170.        (cdr f1)))))
  171.  
  172. (defun f- (a1 &optional a2)        ; unary or binary minus
  173.   "Returns the difference of two floating point numbers."
  174.   (if a2
  175.       (f+ a1 (f- a2))
  176.     (normalize (cons (- (car a1)) (cdr a1)))))
  177.  
  178. (defun f* (a1 a2)            ; multiply in halfword chunks
  179.   "Returns the product of two floating point numbers."
  180.   (let* ((i1 (car (fabs a1)))
  181.      (i2 (car (fabs a2)))
  182.      (sign (not (same-sign a1 a2)))
  183.      (prodlo (+ (hihalf (* (lohalf i1) (lohalf i2)))
  184.             (lohalf (* (hihalf i1) (lohalf i2)))
  185.             (lohalf (* (lohalf i1) (hihalf i2)))))
  186.      (prodhi (+ (* (hihalf i1) (hihalf i2))
  187.             (hihalf (* (hihalf i1) (lohalf i2)))
  188.             (hihalf (* (lohalf i1) (hihalf i2)))
  189.             (hihalf prodlo))))
  190.     (if (> (lohalf prodlo) round-limit)
  191.     (setq prodhi (1+ prodhi)))    ; round off truncated bits
  192.     (normalize
  193.      (cons (if sign (- prodhi) prodhi)
  194.        (+ (cdr (fabs a1)) (cdr (fabs a2)) mantissa-bits)))))
  195.  
  196. (defun f/ (a1 a2)            ; SLOW subtract-and-shift algorithm
  197.   "Returns the quotient of two floating point numbers."
  198.   (if (zerop (car a2))            ; if divide by 0
  199.       (signal 'arith-error (list "attempt to divide by zero" a1 a2))
  200.     (let ((bits (1- maxbit))
  201.       (quotient 0) 
  202.       (dividend (car (fabs a1)))
  203.       (divisor (car (fabs a2)))
  204.       (sign (not (same-sign a1 a2))))
  205.       (while (natnump bits)
  206.     (if (< (- dividend divisor) 0)
  207.         (setq quotient (ash quotient 1))
  208.       (setq quotient (1+ (ash quotient 1))
  209.         dividend (- dividend divisor)))
  210.     (setq dividend (ash dividend 1)
  211.           bits (1- bits)))
  212.       (normalize
  213.        (cons (if sign (- quotient) quotient)
  214.          (- (cdr (fabs a1)) (cdr (fabs a2)) (1- maxbit)))))))
  215.   
  216. (defun f% (a1 a2)
  217.   "Returns the remainder of first floating point number divided by second."
  218.   (f- a1 (f* (ftrunc (f/ a1 a2)) a2)))
  219.       
  220.  
  221. ;; Comparison functions
  222. (defun f= (a1 a2)
  223.   "Returns t if two floating point numbers are equal, nil otherwise."
  224.   (equal a1 a2))
  225.  
  226. (defun f> (a1 a2)
  227.   "Returns t if first floating point number is greater than second,
  228. nil otherwise."
  229.   (cond ((and (natnump (car a1)) (< (car a2) 0)) 
  230.      t)                ; a1 nonnegative, a2 negative
  231.     ((and (> (car a1) 0) (<= (car a2) 0))
  232.      t)                ; a1 positive, a2 nonpositive
  233.     ((and (<= (car a1) 0) (natnump (car a2)))
  234.      nil)                ; a1 nonpos, a2 nonneg
  235.     ((/= (cdr a1) (cdr a2))        ; same signs.  exponents differ
  236.      (> (cdr a1) (cdr a2)))        ; compare the mantissas.
  237.     (t
  238.      (> (car a1) (car a2)))))    ; same exponents.
  239.  
  240. (defun f>= (a1 a2)
  241.   "Returns t if first floating point number is greater than or equal to 
  242. second, nil otherwise."
  243.   (or (f> a1 a2) (f= a1 a2)))
  244.  
  245. (defun f< (a1 a2)
  246.   "Returns t if first floating point number is less than second,
  247. nil otherwise."
  248.   (not (f>= a1 a2)))
  249.  
  250. (defun f<= (a1 a2)
  251.   "Returns t if first floating point number is less than or equal to
  252. second, nil otherwise."
  253.   (not (f> a1 a2)))
  254.  
  255. (defun f/= (a1 a2)
  256.   "Returns t if first floating point number is not equal to second,
  257. nil otherwise."
  258.   (not (f= a1 a2)))
  259.  
  260. (defun fmin (a1 a2)
  261.   "Returns the minimum of two floating point numbers."
  262.   (if (f< a1 a2) a1 a2))
  263.  
  264. (defun fmax (a1 a2)
  265.   "Returns the maximum of two floating point numbers."
  266.   (if (f> a1 a2) a1 a2))
  267.       
  268. (defun fzerop (fnum)
  269.   "Returns t if the floating point number is zero, nil otherwise."
  270.   (= (car fnum) 0))
  271.  
  272. (defun floatp (fnum)
  273.   "Returns t if the arg is a floating point number, nil otherwise."
  274.   (and (consp fnum) (integerp (car fnum)) (integerp (cdr fnum))))
  275.  
  276. ;; Conversion routines
  277. (defun f (int)
  278.   "Convert the integer argument to floating point, like a C cast operator."
  279.   (normalize (cons int '0)))
  280.  
  281. (defun int-to-hex-string (int)
  282.   "Convert the integer argument to a C-style hexadecimal string."
  283.   (let ((shiftval -20)
  284.     (str "0x")
  285.     (hex-chars "0123456789ABCDEF"))
  286.     (while (<= shiftval 0)
  287.       (setq str (concat str (char-to-string 
  288.             (aref hex-chars
  289.                   (logand (lsh int shiftval) 15))))
  290.         shiftval (+ shiftval 4)))
  291.     str))
  292.  
  293. (defun ftrunc (fnum)            ; truncate fractional part
  294.   "Truncate the fractional part of a floating point number."
  295.   (cond ((natnump (cdr fnum))        ; it's all integer, return number as is
  296.      fnum)
  297.     ((<= (cdr fnum) (- maxbit))    ; it's all fractional, return 0
  298.      '(0 . 1))
  299.     (t                ; otherwise mask out fractional bits
  300.      (let ((mant (car fnum)) (exp (cdr fnum)))
  301.        (normalize 
  302.         (cons (if (natnump mant)    ; if negative, use absolute value
  303.               (ash (ash mant exp) (- exp))
  304.             (- (ash (ash (- mant) exp) (- exp))))
  305.           exp))))))
  306.  
  307. (defun fint (fnum)            ; truncate and convert to integer
  308.   "Convert the floating point number to integer, with truncation, 
  309. like a C cast operator."
  310.   (let* ((tf (ftrunc fnum)) (tint (car tf)) (texp (cdr tf)))
  311.     (cond ((>= texp mantissa-bits)    ; too high, return "maxint"
  312.        mantissa-maxval)
  313.       ((<= texp (- mantissa-bits))    ; too low, return "minint"
  314.        mantissa-minval)
  315.       (t                ; in range
  316.        (ash tint texp)))))        ; shift so that exponent is 0
  317.  
  318. (defun float-to-string (fnum &optional sci)
  319.   "Convert the floating point number to a decimal string.
  320. Optional second argument non-nil means use scientific notation."
  321.   (let* ((value (fabs fnum)) (sign (< (car fnum) 0))
  322.      (power 0) (result 0) (str "") 
  323.      (temp 0) (pow10 _f1))
  324.  
  325.     (if (f= fnum _f0)
  326.     "0"
  327.       (if (f>= value _f1)            ; find largest power of 10 <= value
  328.       (progn                ; value >= 1, power is positive
  329.         (while (f<= (setq temp (f* pow10 highest-power-of-10)) value)
  330.           (setq pow10 temp
  331.             power (+ power decimal-digits)))
  332.         (while (f<= (setq temp (f* pow10 _f10)) value)
  333.           (setq pow10 temp
  334.             power (1+ power))))
  335.     (progn                ; value < 1, power is negative
  336.       (while (f> (setq temp (f/ pow10 highest-power-of-10)) value)
  337.         (setq pow10 temp
  338.           power (- power decimal-digits)))
  339.       (while (f> pow10 value)
  340.         (setq pow10 (f/ pow10 _f10)
  341.           power (1- power)))))
  342.                       ; get value in range 100000 to 999999
  343.       (setq value (f* (f/ value pow10) all-decimal-digs-minval)
  344.         result (ftrunc value))
  345.       (let (int)
  346.     (if (f> (f- value result) _f1/2)    ; round up if remainder > 0.5
  347.         (setq int (1+ (fint result)))
  348.       (setq int (fint result)))
  349.     (setq str (int-to-string int))
  350.     (if (>= int 1000000)
  351.         (setq power (1+ power))))
  352.  
  353.       (if sci                ; scientific notation
  354.       (setq str (concat (substring str 0 1) "." (substring str 1)
  355.                 "E" (int-to-string power)))
  356.  
  357.                       ; regular decimal string
  358.     (cond ((>= power (1- decimal-digits))
  359.                       ; large power, append zeroes
  360.            (let ((zeroes (- power decimal-digits)))
  361.          (while (natnump zeroes)
  362.            (setq str (concat str "0")
  363.              zeroes (1- zeroes)))))
  364.  
  365.                       ; negative power, prepend decimal
  366.           ((< power 0)        ; point and zeroes
  367.            (let ((zeroes (- (- power) 2)))
  368.          (while (natnump zeroes)
  369.            (setq str (concat "0" str)
  370.              zeroes (1- zeroes)))
  371.          (setq str (concat "0." str))))
  372.  
  373.           (t                ; in range, insert decimal point
  374.            (setq str (concat
  375.               (substring str 0 (1+ power))
  376.               "."
  377.               (substring str (1+ power)))))))
  378.  
  379.       (if sign                ; if negative, prepend minus sign
  380.       (concat "-" str)
  381.     str))))
  382.  
  383.     
  384. ;; string to float conversion.
  385. ;; accepts scientific notation, but ignores anything after the first two
  386. ;; digits of the exponent.
  387. (defun string-to-float (str)
  388.   "Convert the string to a floating point number.
  389. Accepts a decimal string in scientific notation, 
  390. with exponent preceded by either E or e.
  391. Only the 6 most significant digits of the integer and fractional parts
  392. are used; only the first two digits of the exponent are used.
  393. Negative signs preceding both the decimal number and the exponent
  394. are recognized."
  395.  
  396.   (if (string-match floating-point-regexp str 0)
  397.       (let (power)
  398.     (f*
  399.      ; calculate the mantissa
  400.      (let* ((int-subst (extract-match str 2))
  401.         (fract-subst (extract-match str 4))
  402.         (digit-string (concat int-subst fract-subst))
  403.         (mant-sign (equal (extract-match str 1) "-"))
  404.         (leading-0s 0) (round-up nil))
  405.  
  406.        ; get rid of leading 0's
  407.        (setq power (- (length int-subst) decimal-digits))
  408.        (while (and (< leading-0s (length digit-string))
  409.                (= (aref digit-string leading-0s) ?0))
  410.          (setq leading-0s (1+ leading-0s)))
  411.        (setq power (- power leading-0s)
  412.          digit-string (substring digit-string leading-0s))
  413.        
  414.        ; if more than 6 digits, round off
  415.        (if (> (length digit-string) decimal-digits)
  416.            (setq round-up (>= (aref digit-string decimal-digits) ?5)
  417.              digit-string (substring digit-string 0 decimal-digits))
  418.          (setq power (+ power (- decimal-digits (length digit-string)))))
  419.  
  420.        ; round up and add minus sign, if necessary
  421.        (f (* (+ (string-to-int digit-string)
  422.             (if round-up 1 0))
  423.          (if mant-sign -1 1))))
  424.        
  425.      ; calculate the exponent (power of ten)
  426.      (let* ((expt-subst (extract-match str 9))
  427.         (expt-sign (equal (extract-match str 8) "-"))
  428.         (expt 0) (chunks 0) (tens 0) (exponent _f1)
  429.         (func 'f*))
  430.  
  431.        (setq expt (+ (* (string-to-int
  432.                  (substring expt-subst 0
  433.                     (min expt-digits (length expt-subst))))
  434.                 (if expt-sign -1 1))
  435.              power))
  436.        (if (< expt 0)        ; if power of 10 negative
  437.            (setq expt (- expt)    ; take abs val of exponent
  438.              func 'f/))        ; and set up to divide, not multiply
  439.  
  440.        (setq chunks (/ expt decimal-digits)
  441.          tens (% expt decimal-digits))
  442.        ; divide or multiply by "chunks" of 10**6
  443.        (while (> chunks 0)    
  444.          (setq exponent (funcall func exponent highest-power-of-10)
  445.            chunks (1- chunks)))
  446.        ; divide or multiply by remaining power of ten
  447.        (funcall func exponent (aref powers-of-10 tens)))))
  448.           
  449.     _f0))                ; if invalid, return 0
  450.  
  451.  
  452.